Robustness of a Network of Networks 
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Abstract 

Almost all network research has been focused on the properties of a single network that does not interact and depends on 
other networks. In reality, many real-world networks interact with other networks. Here we develop an analytical framework for 
studying interacting networks and present an exact percolation law for a network of n interdependent networks. In particular, we 
find that for n Erdos-Renyi networks each of average degree k, the giant component, Poo, is given by Poo — p[l — exp(— fcPoo)]™ 
where 1 — p is the initial fraction of removed nodes. Our general result coincides for n = 1 with the known Erdos-Renyi 
second-order phase transition for a single network. For any n > 2 cascading failures occur and the transition becomes a 
first-order percolation transition. The new law for Pec shows that percolation theory that is extensively studied in physics and 
mathematics is a limiting case (n = 1) of a more general general and different percolation law for interdependent networks. 



In recent years dramatic advances in the field of com- 
plex networks have occurred (ll-fl3j. The internet, airline 
routes, and electric power grids are all examples of net- 
works whose function relies crucially on the connectivity 
between the network components. An important prop- 
erty of such systems is their robustness to node failures. 
Almost all research has been concentrated on the case 
of a single or isolated network which does not interact 
with other networks. Recently, based on the motivation 
that modern infrastructures are becoming significantly 
more dependent on each other, a system of two coupled 
interdependent networks has been studied [lH • A funda- 
mental property of interdependent networks is that when 
nodes in one network fail, they may lead to the failure 
of dependent nodes in other networks which may cause 
further damage in the first network and so on, leading to 
a global cascade of failures. Buldyrev et al. [l4| devel- 
oped a framework for analyzing robustness of two inter- 
acting networks subject to such cascading failures. They 
found that interdependent networks become significantly 
more vulnerable compared to their noninteracting coun- 
terparts. 

For many important examples, more than two net- 
works depend on each other. For example, diverse infras- 
tructures are coupled together, such as water and food 
supply, communications, fuel, financial transactions, and 
power stations [l5l - [l8j . For further examples see Section 
I in the Supplementary Information (SI). Understanding 
the robustness due to such interdependencies is one of the 
major challenges for designing resilient infrastructures. 

We introduce here a model system, comprising a net- 
work of n coupled networks, where each network con- 
sists of N nodes (See Fig. 1). The N nodes in each 
network are connected to nodes in neighboring networks 
by bidirectional dependency links, thereby establishing a 
one-to-one correspondence as illustrated in Fig. 2 in SI. 
We develop a mathematical framework to study the ro- 



bustness of this "network of networks" (NON). We find 
an exact analytical law for percolation of a NON system 
composed of n coupled randomly connected networks. 
Our result generalizes the known Erdos-Renyi (ER) frol - 
ic result for the giant component of a single network 
and the n = 2 result found recently [lij, and shows that 
while for n = 1 the percolation transition is a second or- 
der transition, for n > 1 cascading failures occur and the 
transition becomes a first order transition. Our results 
for n interdependent networks suggest that the classi- 
cal percolation theory extensively studied in physics and 
mathematics is a limiting case of a general theory of per- 
colation in NONs, or networks with multiple types of 
connectivity links. This general theory has many novel 
features which are not present in classical percolation. 
Additionally, we find: 

(i) the robustness of NON significantly decreases with 
n, and 

(ii) for a network of n ER networks all with the same 
average degree k, there exists a minimum degree fc m i n (n) 
increasing with n, below which p c = 1, i.e., for k < fc m i n 
the NON will collapse once any finite number of nodes 
fail. We find an analytical expression for k m i n (n), which 
generalizes the known result fc m in(l) = 1 for ER below 
which the network collapses. We also discuss the critical 
effect of loops in the NON structure. 

Real-world interacting networks (See SI for more de- 
tails) arc characterized by complex correlations and a va- 
riety of organizational principles governing their internal 
structures and interdependencies. Once these correla- 
tions are quantified from the statistical analysis of actual 
data bases and the organizational principles are specified 
from the engineering literature, real world networks can 
be studied by computer simulations. These simulations 
will have many parameters and therefore their outcome 
will also require complex interpretation. It is therefore 
very important to develop simple analytically tractable 
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models for the robustness of interdependent networks 
against which such simulations can be tested. Well- 
known examples of simplified models that both demon- 
strate a fundamental phenomenon and significantly ad- 
vance our knowledge are the Ising model in statistical 
mechanics and the Erdos-Renyi model in graph theory. 
This paper presents a simple model that can serve as a 
benchmark for further studies of NONs. 

We assume that a network m (m = 1,2, ...,n) in the 
NON is a randomly connected network with a degree dis- 
tribution P m (k). We call a pair of networks A and B a 
fully interdependent pair if it satisfies the following con- 
dition: each node Ai of network A is connected to one 
and only one node Bi in network B by a bidirectional de- 
pendency link such that if node Ai fails, Bi also fails and 
vice versa. Since the number of nodes in each network 
is the same, these bidirectional links establish a one-to- 
one correspondence between the nodes in the networks 
belonging to an interdependent pair. Each node of the 
NON represents a network and each edge represents a 
fully interdependent pair of networks. First, we will dis- 
cuss the case when the NON is a loopless tree of n net- 
works (Fig. 1). The dependency edges in such a NON es- 
tablish a unique one-to-one correspondence between the 
nodes of any two networks not necessarily belonging to 
the same fully interdependent pair. This one-to-one cor- 
respondence established by the intcrdependency links be- 
tween the nodes of different networks in the loopless NON 
uniquely maps any set of nodes in one of the networks to 
a set of nodes (which we will call an image of the original 
set) in any other network of the NON (See SI for more 
details). In principle, the assumption of full interdepen- 
dence allows one to collapse all the networks of the NON 
onto a single network with multiple types of links. 

We assume that in order to remain functional a node 
must belong to a sufficiently large mutually connected 
cluster [l4| (See detailed definition in SI). We will show 
that a large mutually connected cluster which includes a 
finite fraction of the nodes in each network exists only in 
networks of sufficiently high mean degree. We call this 
large mutually connected cluster a mutual giant compo- 
nent, and we postulate that only nodes in the mutual 
giant component remain functional. 

We assume that due to an attack or random failure 
only a fraction of nodes p in one particular network which 
we will call the root of the NON. We can now observe a 
cascade of failures caused by the failure of the dependent 
nodes in the networks connected directly to the root by 
the edges of the NON. The damage will further spread 
to more distant networks. Moreover, fragmentation of 
each network caused by the removal of certain nodes will 
cause malfunction of other nodes which will now belong 
to small isolated clusters. This malfunction will cause 
dependent nodes in neighboring networks to malfunction 
as well. Depending on the time scales of these processes, 
the damage can spread across the NON back and forth, 
which we can visualize as cascades of failures, as shown 
in Fig. 3 of the SI section. At the end of this process only 



the mutual giant component of the NON, if it exists, will 
remain functional. 

We now introduce generating functions [3, [22| - [25j 
of each network, G m o{z) = ^ P m {k)z k , and the gen- 
erating function of the associated branc hing process, 
G m i(z) = G' m0 {z)/G' m0 (l). It is known Hi [23[ that 
the generating functions of a randomly connected net- 
work once a fraction 1 — p of nodes has been ran- 
domly removed are the same as the generating func- 
tions of the original network with the new argument 
1 — p(l — z). Furthermore it is known [24|, [25[ that 
the fraction of nodes in the giant component of a single 
randomly connected network is = pg m {p), where 

g m {p) = 1 - G m0 (l -p(l - f m )) and f m satisfies a tran- 
scendental equation f m (p) = G m i(l - p(l - f m {p)))- 

We next prove that the fraction of nodes, /ioo,n, in 
the mutual giant component of a NON composed of n 
networks is the product: 

n 

Moo.n = P Y\_ 9m(Xm), (1) 
m— 1 

where each x m satisfies the equation 

— l^oo^n/Qm (^-m ) ■ (2) 

The system of n + 1 equations (fT]) and @ defines n + 1 
unknowns: fx OCjn ,xi,X2, ■■■,x n as functions of p and the 
degree distributions P m {k). 

We derive Eqs. (1) and (2) by mathematical induc- 
tion. (An alternative proof is given in the SI) . Indeed, for 
n = 1, Eqs. (1) and (2) follow directly from the definition 
of Hoo.i ■ Assuming that the formulas are valid for a NON 
of n — 1 networks we will prove that they are valid also for 
a NON of n networks. A loopless NON of n networks can 
be always represented as one of its networks connected by 
a single edge to the other n — 1 networks in the NON. All 
the nodes in the 77-th network, which do not belong to the 
image of the mutual giant component yUoo,n-i of the NON 
of 77 — 1 networks will stop to function. The fraction of the 
nodes in the image of this mutual giant component onto 
the n-th network satisfies the equation x\, n — /ioo.n-i(p)- 
The fraction of nodes belonging to the giant component 
of this dependency image is fii t7l = xi, n g n {xi tn ). Only 
the nodes in the NON of n — 1 networks which belong 
to the dependency image of the giant component of the 
n-th network will remain functional. Due to the random- 
ness of the connectivity links in different networks, this 
dependency image can be represented as a random se- 
lection of the fraction g n {xi tn ) out of the originally sur- 
vived nodes, or as random selection of p\ = pg n (xi,n) 
fraction of nodes in one of the networks comprising the 
NON of n — 1 networks. The fraction of nodes in the 
new mutual giant component of the NON of 77 — 1 net- 
works corresponding to this new random selection will be 
i(pi). The image of this mutual giant component 
in the 77-th network is equivalent to a random selection of 
x 2,n = ^oc.n~i{pi) / 9n{xi. n ) fraction of nodes out of the 
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entire n-th network. As we continue this process, the se- 
quence of giant components fj,j in in the n-th network, ran- 
domly selected sets x n in the n-th network and randomly 
selected sets pj in the NON of n — 1 networks will satisfy 
the recursion relations Xj+i >n = (J*oo,n—i(pj)/9n(%j,n)i 

(J-j+l,n = Xj + i in g n (Xj + l t n), Pj+i = pg n (Xj + i tn ). 

In the limit j — > oo, this process will converge, i.e. all 
the parameters in the two successive steps will coincide: 

x j+l,n ^ Xj,n = Xn, Pj ^ P9n(Xn) and /ioo,n— 1 

(Pi) "> 

fJ-<x>,n- Then x n = fJ-oo,n/ g-n(x n ) which is identical 
to the last equation in Eqs. (2) and fJioo,n-i(Pj) ~ 5 " 

P9n(x n )Y[^n=l9m.(x m ) = P Tim=l 9m{x m ) = which 

is identical to Eq. (1). By the assumption of induction 
x m = fJ,oo,n-i(pj) 1 9m(xm) = Moo, n/9m(x m ) which com- 
pletes the set of Eqs. (2). Finally fij+i, n — > x n g n (x n ) = 
Hoo.n/ 9n{x n ) and thus the fraction of nodes in the giant 
n-th network coincides with the mutual giant component 
in the NON of n — 1 networks. The SI presents an al- 
ternative analytical derivation of Eqs. (1) and (2), which 
represent a certain type of cascading failures. The SI also 
presents simulation results which agree well with the the- 
ory (Figs. 5 and 6 in SI). 

For the case of a NON with loops, the closed path of 
fully interdependent pairs starting form a network A and 
ending at the same network A will establish a dependence 
of nodes Ai on node Aj t , where ji is a transposition of i. 
Then the failure of single node i will cause an entire cycle 
in the transposition to fail. The average size of a cycle 
in the transposition of N elements grows as N / In AT, so 
the initial failure of In N nodes will cause almost all the 
nodes of the NON to fail without taking into account 
any connectivity links which will cause additional dam- 
age. So the NON with loops is unstable against removal 
of an infinitely small fraction of nodes unless the transpo- 
sition ji is not random. In case when the transposition 
ji is trivial, ji = i, we have the same one-to-one cor- 
respondence between the nodes as in the loopless NON 
and then Eq. ([T]) and @ arc valid. This is since in our 
proof we did not use any other property of a NON ex- 
cept the unique one-to-one correspondence of the nodes 
in different networks. 

The case of NON of n Erdos-Renyi (ER) (HHI 
networks with average degrees k\, k%, ...fcj, k„ can be 
solved explicitly. In this case, we have G\ t i(x) = 
Go,i(x) = exp[ki(x - 1)] H3]. Accordingly gi{x t ) = 
1 - cxp[fciXi(/i - 1)], where fc = exjp[kiXi(fc - 1)] and 
thus gi(xi) = 1 — fc. Using Eq. ([2]) for Xi we get 

n 

fc = cxp[-pki JJ(1 - fc)],i - 1,2,..., n. (3) 

j=i 

These equations can be solved analytically, as shown in 
detail in the SI section. They have only a trivial solution 
{fc = 1) if p < p c , where p c is the mutual percolation 
threshold. When the n networks have the same average 
degree fc, ki — k (i = 1, 2, ...,n), we obtain from Eq. (3) 
that f c = fc(p c ) satisfies 

/c = ew. (4) 



where the solution can be expressed in term of the Lam- 
bert function W(x) HHEl, fc = -[nW(-^e-i)]- 1 . 

Once fc is known, we obtain p c and //oo.n = -Poo by 
substituting fc,; = k into Eq. (S10) of the SI section 

Pc = [nkf c {l-fc) {n - 1) ]-\ 

P _ We • ( 5 ) 

r °° ~ nkf a ■ 

For n = 1 we obtain the known result p c = 1/fc of Erdos- 
Renyi [l9l42l| . Substituting n = 2 in Eqs. dU) and © 
one obtains the exact results of [Til ]. 

To analyze p c as a function of n, we find f c from Eq. 
(|4]) and substitute it into Eq. ([5]), and we obtain p c as a 
function of n for different k values, as shown in Fig. 2(a). 
It is seen that the NON becomes more vulnerable with in- 
creasing n or decreasing fc (p c increases when n increases 
or fc decreases). Furthermore, for a fixed n, when fc is 
smaller than a critical number fc m j n (n), p c > 1 meaning 
that for fc < k m i n (n), the NON will collapse even if a 
single node fails. Fig. 2 (b) shows the minimum average 
degree fc m ; n as a function of the number of networks n. 
From Eq. ([S]) we get the minimum of fc as a function of 
n 

femin(n) = [n/ c (l - /c)^]" 1 , (6) 

Note that Eq. together with Eq. Q yield the value 
of fcmin(l) = 1, reproducing the known ER result, that 
(fc) = 1 is the minimum average degree needed to have 
a giant com pon ent. For n = 2, Eq. ^ yields the result 
obtained in [1J], i.e., fc m i n = 2.4554. 

From Eqs. (l)-(3) we obtain the percolation law for the 
order parameter, the size of the mutual giant component 
for all p values and for all fc and n, 

Moo.n = Poo = P[l - CXp(-fcP 00 )]™. (7) 

The solutions of equation (JT]) arc shown in Fig. 3 for 
several values of n. Results are in excellent agreement 
with simulations. The special case n = 1 is the known 
ER percolation law for a single network [19I - I2TI ]. Note 
that Eqs. (4)-(7) are based on the assumption that all n 
networks have the same average degree fc. 

In summary, we have developed a framework, Eqs. (1) 
and (2), for studying percolation of NON from which 
we derived an exact analytical law, Eq. ([?]). for percola- 
tion in the case of a network of n coupled ER networks. 
Equation {?]) represents a bound for the case of partially 
dependent networks [28| . which will be more robust. In 
particular for any n > 2, cascades of failures naturally ap- 
pear and the phase transition becomes first order transi- 
tion compared to a second order transition in the classical 
percolation of a single network. These findings show that 
the percolation theory of a single network is a limiting 
case of a more general case of percolation of interdepen- 
dent networks. Due to cascading failures which increase 
with n, vulnerability significantly increases with n. We 
also find that for any loopless network of networks the 
critical percolation threshold and the mutual giant com- 
ponent depend only on the number of networks and not 
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on the topology (see Fig. 1(a)). When the NON includes mutual giant component exists, 
loops, and dependency links are random, p c = 1 and no 



[1] Watts D. J. & Strogatz S. H. Nature 393, 440-442 (1998). 
[2] Albert R., Jeong H. & Barabasi A. L. Nature 406, 378- 
382 (2000). 

[3] Cohen R. et al. Phys. Rev. Lett. 85, 4626-4628 (2000). 
[4] Callaway D. S. et al. Phys. Rev. Lett. 85, 5468-5471 
(2000). 

[5] Albert R. & Barabasi A. L. Rev. Mod. Phys. 74, 47-97 
(2002). 

[6] Newman M. E. J. SIAM Review 45, 167-256 (2003). 

[7] Dorogovtsev S. N. & Mendes J. F. F. Evolution of Net- 
works: From Biological Nets to the Internet and WWW 
(Physics) (Oxford Univ. Press, New York, 2003). 

[8] Song C. et al. Nature 433, 392-395 (2005). 

[9] Satorras R. P. & Vespignani A. Evolution and Structure 
of the Internet: A Statistical Physics Approach (Cam- 
bridge Univ. Press, England, 2006). 
[10] Caldarelli G. & Vespignani A. Large scale Structure and 

Dynamics of Complex Webs (World Scientific, 2007). 
[11] Barrat A., Barthelemy M. & Vespignani A. Dynamical 
Processes on Complex Networks (Cambridge Univ. Press, 
England, 2008). 
[12] Havlin S. & Cohen R. Complex Networks: Structure, Ro- 
bustness and Function (Cambridge Univ. Press, England, 
2010). 

[13] Newman M. E. J. Networks: An Introduction (Oxford 



Univ. Press, New York, 2010). 
[14] Buldyrev S. V. et al. Nature 464, 1025-1028 (2010). 
[15] Peerenboom J., Fischer R. & Whitfield R. in Pro. 

CRIS/DRM/IIIT/NSF Workshop Mitigat. Vulnerab. 

Crit. Infrastruct. Catastr. Failures (2001). 
[16] Rinaldi S., Peerenboom J. & Kelly T. IEEE Contr. Syst. 

Mag. 21, 11-25 (2001). 
[17] Rosato V. et al. Int. J. Crit. Infrastruct. 4, 63-79 (2008). 
[18] Vespignani A. Nature 464, 984-985 (2010). 
[19] Erdos P. & Renyi A. I. Publ. Math. 6, 290-297 (1959). 
[20] Erdos P. & Renyi A. Publ. Math. Inst. Hung. Acad. Sci. 

5, 17-61 (1960). 
[21] Bollobas B. Random Graphs (Academic, London, 1985). 
[22] Newman M. E. J. Strogatz S. H. & Watts D. J., Phys. 

Rev. E 64, 026118 (2001). 
[23] Newman M. E. J. Phys. Rev. E 66, 016128 (2002). 
[24] Shao J. et al. Europhys. Lett. 84, 48004 (2008). 
[25] Shao J. et al. Phys. Rev. E 80, 036105 (2009). 
[26] Lambert J. H. Acta Helveticae physico mathematico 

anatomico botanico medica, Band III, 128-168, (1758). 
[27] Corless R. M. et al. Adv. Computational Maths. 5, 329- 

359 (1996). 

[28] Parshani R. et al. Phys. Rev. Lett. 105, 048701 (2010). 



4 




line-like NON star-like NON tree-like NON 



FIG. 1: (color online) Three types of loopless networks 
of networks composed of five coupled networks. All have 
same percolation threshold and same giant component. 
The darker green node is the origin network on which 
failures occur. 
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FIG. 2: (a) The critical fraction p c for different k and n 
and (b) minimum average degree k mm as a function of 
the number of networks n. The results of (a) and (b) are 
obtained using Eqs. ([5]) and (|6]) respectively and are in 
good agreement with simulations. In simulations p c was 
calculated from the number of cascading failures which 
diverge at p c [28[ (see also Fig. 7 in SI). 
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FIG. 3: For looplcss NON, as a function of p for 
k = 5 and several values of n. The results obtained using 
Eq. ([7]) agree well with simulations. 
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